#version 430 core

#define LOC_X 16
#define LOC_Y 16

#define GRID_SIZE_X 64
#define GRID_SIZE_Y 64

layout(local_size_x=LOC_X, local_size_y=LOC_Y) in;

layout(binding=0) uniform sampler2D levelSetImage;
layout(binding=1, r32i) uniform writeonly iimage2D nodeIndexImage;

layout(binding=0, offset=4) uniform atomic_uint counter;

layout(std430, binding=0) coherent buffer femElementPositions {
	ivec2 elementPositions[];
};

void main(void)
{
	// one invocation per cell
	ivec2 gid = ivec2(gl_GlobalInvocationID.xy);
	bvec2 cond = lessThan(gid, ivec2(GRID_SIZE_X, GRID_SIZE_Y));
	if (cond.x && cond.y)
	{
		ivec2 off1 = gid + ivec2(1,0);
		ivec2 off2 = gid + ivec2(0,1);
		ivec2 off3 = gid + ivec2(1,1);

		// level set values at the four corners
		vec4 levelSetValues = textureGather(levelSetImage, (gl_GlobalInvocationID.xy+vec2(1.0,1.0))*vec2(1.0/(GRID_SIZE_X+1.0), 1.0/(GRID_SIZE_Y+1.0)));

		ivec4 signLevelSetValues = ivec4(sign(levelSetValues));
		int sum = signLevelSetValues.x + signLevelSetValues.y + signLevelSetValues.z + signLevelSetValues.w;
	
		if (sum != 4)
		{
			// cell is not fully air and is part of the FEM calculation

			// mark the four corner nodes with index 0 (other nodes will remain -1)
			imageStore(nodeIndexImage, gid , ivec4(0,0,0,1));
			imageStore(nodeIndexImage, off1, ivec4(0,0,0,1));
			imageStore(nodeIndexImage, off2, ivec4(0,0,0,1));
			imageStore(nodeIndexImage, off3, ivec4(0,0,0,1));
			
			// FEM cell is assigned a unique number, save integer coordinates in buffer
			uint femIndex = atomicCounterIncrement(counter);
			elementPositions[femIndex] = gid;
		}
	}
}